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Collapse and reverse to collapse explosion transition in self-gravitating systems are studied by 
molecular dynamics simulations. A microcanonical ensemble of point particles confined to a spherical 
box is considered; the particles interact via an attractive soft Coulomb potential. It is observed that 
the collapse in the particle system indeed takes place when the energy of the uniform state is put 
near or below the metastability-instability threshold (collapse energy), predicted by the mean- field 
theory. Similarly, the explosion in the particle system occurs when the energy of the core-halo 
state is increased above the explosion energy, where according to the mean field predictions the 
core-halo state becomes unstable. For a system consisting of f 25 ~ 500 particles, the collapse takes 
about fO'' single particle crossing times to complete, while a typical explosion is by an order of 
magnitude faster. A finite lifetime of metastable states is observed. It is also found that the mean- 
field description of the uniform and the core-halo states is exact within the statistical uncertainty 
of the molecular dynamics data. 

PACS numbers: 64.60.-i 02.30.Rz 04.40.-b 05.70.Fh 



INTRODUCTION 

Systems of particles interacting via a potential with at- 
tractive nonintegrable large r asymptotics, U{r) ^ r~", 
< a < 3, and sufficiently short-range small r regu- 
larization exhibit gravitational phase transition between 
a relatively uniform high energy state and a low-energy 
state with a core-halo structure [l S 0, Sll 0,11 113 ■ Ex- 
tensive mean-field (MF) studies of the equilibrium prop- 
erties of such systems [100|lll011Il3| revealed 
that in a microcanonical ensemble during such a transi- 
tion entropy has to undergo a discontinuous jump from 
a state that just ceases to be a local entropy maximum 
to a state with the same energy but different tempera- 
ture, which is the global entropy maximum. Due to the 
long-range nature of gravitational interaction, the MF 
studies are believed to provide asymptotically (in the in- 
finite system limit) exact information about the density 
and the velocity distributions and other thermodynami- 
cal parameters of the uniform state. The applicability of 
the MF theory to the description of the core-halo state 
is less obvious as the properties of a core are controlled 
by the short-range asymptotics of the potential. 

Relatively little is known about how such a transi- 
tion actually occurs, however. Youngkins and Miller 4] 
performed a Molecular Dynamics (MD) study of a one- 
dimensional system consisting of concentric spherical 
shells. Their main emphasis was to check the MF de- 
scription of the stable and metastable states rather than 
to study the dynamics of the phase transition itself. 
Cerruti-Sola, Cipriani, and Pettini studied the phase 
diagram of a more realistic 3-dimensional particle system 
by using Monte Carlo and MD methods. Their studies 
again focused on the equilibrium properties of the system 
rather than on the dynamics of the transitions. In addi- 



tion, their general conclusions about the second order of 
the gravitational phase transition apparently contradict 
the MF results |lli,|i,[l3- 

Here, we attempt to resolve this contradiction. A MF 
description of the dynamics of collapse in ensembles of 
self-gravitating Brownian particles with a bare 1/r in- 
teraction based on a Smoluchowski equation was devel- 
oped by Chavanis et al It predicts a self-similar 
evolution of the central part of density distribution to 
a finite-time singularity. However, the precise nature of 
the random force and friction terms in the corresponding 
Fokkcr-Plank equation as well as the applicability of the 
overdamped limit used to reduce the Fokker-Plank equa- 
tion to the Smoluchowski equation are not entirely under- 
stood. A more rigorous approach based on the Fokker- 
Plank equation with Landau collision integral was used 
by Lancellotti and Kiessling to prove a scaling prop- 
erty of the central part of the density profile. The model 
considered there allows the particles to escape to infin- 
ity and therefore does not have an equilibrium or even a 
metastable state. 

There exists a vast amount of literature on 
cosmologically- and astrophysically- motivated studies 
of the temporal evolution of naturally occurring self- 
gravitating systems (see, e.g., Ref. 14] and references 
therein). The selection of systems and their initial 
and final conditions made in such studies are typically 
astrophysically-motivated; the considered systems are of- 
ten too complex to make a general conclusions about the 
phase diagrams and phase transitions in such systems. 

In this paper we present MD studies of gravitational 
collapse, and reverse to collapse, i.e., explosion, transi- 
tions in a microcanonical ensemble of self-attracting par- 
ticles. Besides their pure statistical mechanical impli- 
cations, these studies represent our attempt to bridge a 
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FIG. 1: Plots of entropy s(e) (solid line) and inverse tem- 
perature /3(e) = ds/de (dashed line) vs. energy e for a system 
with a gravitational phase transition and a short-range cutoff. 



gap between the usually complicated MD and hydrody- 
namic simulations of the realistic astrophysical systems 
and the MF analysis of the phase diagram of simple self- 
gravitating models. 



A system with soft Coulomb potential —{P + ?'o) 



-1/2 

is considered. Such systems are have been studied us- 
ing both MF theory (see, e.g., Refs. 0,01) and simu- 
lations We chose the microcanonical ensemble as 
the most fundamental one for the long-range interacting 
systems. It has to be noted that the considered system 
is strongly ensemble-dependent: While the nature of the 
uniform state is the same in both microcanonical and 
canonical ensembles (apart form the difference in their 
stability range), the core- halo states and the collapse it- 
self in these ensembles have very little in common with 
each other 

A MF phase diagram of the considered self-attracting 
microcanonical system is presented in Fig. (Q) 0, lioj . 
High- and low-energy branches terminating at the ener- 
gies tcoii and texpi correspond to the uniform and core- 
halo states. The energy e* where the entropies of the 
core-halo and uniform states are equal is the energy of 
the true phase transition; the uniform and the core-halo 
states are metastable in the energy intervals {tcoiii e*) and 
(e*, ^expi), respectively. However, for the phase transition 
to occur at or near e*, a macroscopic-scale fluctuation 
with prohibitively low entropy is required. Consequently, 
the metastable branches are stable everywhere except at 
ncinitv 

Mm, 

system. Hence it is natural to assume that once the en- 
ergy of the system in the uniform state is set sufficiently 
near above or below ecoih the system will undergo a col- 



the vicinity of Ae ^/'^ of their end-points ecoii and 

^expi [hI , where N is the number of particles in the 



lapse to a core-halo state with the same energy and higher 
entropy. Similarly, if the energy of the core-halo system is 
set sufficiently near below or above eexpi, the system will 
undergo an explosion bringing it to a uniform state with 
the same energy and higher entropy. Our goal here is to 
study if and how such collapse and explosion proceed in a 
realistic three-dimensional N-particle dynamical system. 

The paper is organized as follows. In the next sec- 
tion we formally introduce the system, outline the MF 
analysis and describe the MD setup. Then, we present 
the simulation results for the equilibrium uniform and 
core-halo states and compare them to the MF predic- 
tions. After that we describe and interpret the observed 
dynamics of the collapse and the explosion transitions. 
A discussion of the obtained results concludes the paper. 



SIMULATION 

We consider a system consisting of A^ identical particles 
of unit mass confined to a spherical container of radius 
R with reflecting walls. The particles interact via the 
attractive soft Coulomb pair potential — (r^ -I- rg)~^/^. 
Using a traditional convention for self-gravitating sys- 
tems in which the equilibrium properties of such systems 
become universal, we define rescaled energy e, inverse 
temperature /3, distance x, velocity u, and time r as 




(1) 



The unit of time, often referred to as crossing time, 
[i] — ^1/2 ; is obtained by dividing the unit of length 
R by the unit of velocity ^JN/R. This unit of time is 
also proportional to the period of plasma oscillations in 
a medium with the charge concentration N/ R^ . As this 
time unit has purely kinematic origin, we do not expect 
the evolution of systems having different A^ and R to 
be universal in time r. The evolution, assuming that 
it is collisional, is expected to be universal in the relax- 
ation time Tr = T^^j^ where the factor N/\nN is 
proportional to the number of crossings a typical particle 
needs to change its velocity by a factor of 2 through weak 
Coulomb scattering events. 

The soft core radius xq = r^/R — 5xl0~'^is chosen to 
be well below the critical value Xgr ~ 0.021, above which 
the collapse-explosion transition is replaced by a normal 
first-order phase transition [lO| . 
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FIG. 2: MF density profiles pv.(x) of a uniform state (dashed 
line) and pc.-h{x) of a core-halo state (solid line) for e = ecoii 



The MF theory of the system is described in de- 
tail in Ref . jH . The equilibrium velocity distribution is 
Maxwellian and isothermal, while the equilibrium (sad- 
dle point) density profile p{x) corresponding a stable 
or a metastable state, is a spherically symmetric so- 
lution of the integral equation Q. This equation re- 
places the Poisson-Boltzmann differential equation for 
the self-consistent potential (see, for example, Ref. Q) 
since the interparticle interaction considered here is not 
pure Coulombic. 



F[p(.),x] = exp 



p(x) = pQF[p{.),yi] 
P(x') 



^(x-x')2 + a;2 



1 

^+2 



p(xi)p(x2) 



=d^xid^X2 



V(xi^X2)2+lEf 



(2) 



The equilibrium density profile p{x) obtained from this 
equation is then used to calculate the entropy and the 
pressure. 

The MF phase diagram of the system is presented in 
Fig. 2] The collapse and explosion energies are Ccoii ~ 
—0.339 and texpi ~ 0.267. Examples of the uniform and 
the core-halo density profiles for e — ecoU are shown in 
Fig.H 

In the MD simulations we consider a system consist- 
ing of iV = 125 ~ 500 particles in a spherical container 
of the radius R = 1. All interparticle forces are cal- 
culated directly at each time step dt. This is done to 
avoid any mean-field-like effects inevitably present in any 
truncated multipole or Fourier potential expansion. The 
particle velocities and coordinates are updated according 
to the velocity- Verlet algorithm which is symplectic and 



reversible. 

The system is initialized by randomly distributing par- 
ticles according to a spherically symmetric density pro- 
file; typically the appropriate MF density profile p{x) 
was used. The potential energy ([/) of the initial con- 
figuration is calculated, and the target kinetic energy 
Ek = E ~ U \s determined. The particle velocities are 
randomly generated from some (usually Maxwell) dis- 
tribution with the appropriate square average. Finally, 
the deviation of the total energy from its target value, 
caused by a stochasticity in velocity assignment, is de- 
termined, and the velocities are rescaled to fine-tune the 
total energy. Due to the isotropicity of the random veloc- 
ity assignment, we have always obtained the states with 
sufficiently low total angular momentum which collapsed 
to single-core states rather than to binaries [l^ . 

To implement the reflective boundary condition, at 
each time step the normal components v±^ of the veloci- 
ties of all particles which had escaped from the container 
were reversed. Values of the normal components were 
stored to evaluate the pressure on the wall P. 



2-kRH" 



(3) 



During each simulation run we measured such charac- 
teristics as the kinetic energy e^in = the virial vari- 
able a (dimension of energy) quantifying deviations from 



the virial theorem, a 



_ 3PVR 
?3 , 



(where P is the 



pressure on the wall, V — 47ri?''/3 is volume of the con- 
tainer, and the factor N'^/R rescalcs the volume-pressure 
term to the unit of energy introduced in (^), ratio of 



the velocity moments, 



(which should be 1 for the 



Maxwell distribution), and the number of particles in the 
core Nc of the prescribed radius Xc- For the last measure- 
ment we count the number of particles Ni that are within 
Xc from the ith particle and find the particle which has 
the largest Ni. 

In addition, we measured the histograms of the veloc- 
ity distribution and radial distribution functions, W{u) 
and C{x), respectively. The latter was defined as the 
number of particles in the spherical layer of the radius x 
around each particle, normalized by the volume of such 
layer, disregarding the nonuniformity of the system and 
the boundary effect. 

The measurements of the "scalar" quantities such as 
energy, kinetic energy, pressure, and velocity distribu- 



tion moments were taken in time intervals t„ 



which 



were selected sufficiently long to avoid measuring the un- 
changed configuration repeatedly and sufficiently short 
not to miss the important details of the system evolution. 
We usually pick T,„eas of the order of the uniform density 
sphere crossing time r^ross — which is a half period of 
the oscillation of a particle released with zero velocity at 
the container wall. The histogram data, such as the ve- 
locity distribution and the radial distribution functions, 
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was incremented at each Tmeas and accumulated over a 
longer time period Thtst, Thist ~ 10 - 10^ x r^eas- 

Our attempts to resolve the high-density part of the 
radial density profile of the system turned out to be fruit- 
less due to the strong fluctuations in the positions of this 
part. This fluctuations result in smearing the central 
peak in both core- halo and low energy uniform states. 
Considering the center of mass system of reference does 
not resolve this difficulty, as, despite being dense, the 
core typically contains only 10 - 20% of the total system 
mass (see below) and the positions of the core and the 
center of mass of the system do not usually coincide. 

To control the quality of the simulation, we monitored 
the total energy e and the total angular momentum L. 
We selected timestep dt small enough to keep the total 
energy variation within 0.01% of its initial value, usually 
we used dt = 10~^, or in rescaled units, dr ~ 10^^. 
For such time steps, the relative deviation of the angular 
momentum was within 10~^^. 

All the measurements below are presented in the 
rescaled dimensionless units as defined in Eq. (Q). 
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FIG. 3; Plots of time dependence of kinetic energy e^in (solid 
line) , virial variable a (dashed line) , and total energy e (dotted 
line) for a uniform system of iV = 250 particles at e = —0.3. 



UNIFORM AND CORE-HALO EQUILIBRIUM 
STATES: COMPARISON TO THE MF 



To check our simulation procedure and possibly re- 
solve the apparent contradiction between the MF and 
the particle simulation results 13], we first considered 
the system in what we expected to be a stable or a 
metastable states far away from a transition point. Since 
we were interested in the equilibrium properties, we were 
initiating the MD systems according the corresponding 
MF predictions. It meant that the density profiles were 
seeded according to the MF profiles and the velocities 
were assigned according to the Maxwell distribution. We 
observed that the MF density initiation virtually elimi- 
nates the transitory period, while the method of velocity 
assignment was practically unimportant, provided that 
it gave the correct value for the total kinetic energy. 
For example, it takes a system initialized with a flat 
VF(u) — const about t ~ to evolve to the Maxwell 
distribution. 

A typical plot of the steady state time dependence of 
the kinetic energy, virial variable, and the total energy is 
presented in Fig. |3| 

The comparison between the MD measurements and 
the MF results for the uniform and core-halo states is 
presented in Tables ^ and ^ and reveals a perfect agree- 
ment between these two sets of data. To obtain the 
expression for the MF virial variable 



(yMF = e + efcm(l - 87rp(l)/3) 



(4) 



we write for the pressure at the container wall P = 
2p{x = l)efei„/3 implying that the system is isothermal. 
Since the interparticle potential is not pure Coulombic, 



TABLE I: Equilibrium MD and MF results for a uniform state 
for e = -0.3, iV = 250, and < r < 5000 





MD 


MF 


e 


-0.3 ±5 X lO"'^ 


-0.3 




0.66 ±0.05 


0.644 


O 


± 0.03 


0.012 


19{i;2)2/5{u*) 


1.01 ±0.04 


1 



the virial variable is non-zero. The difference is especially 
prominent for the core-halo states where more particles 
"probe" the short-range part of the potential. 

To evaluate the core radius and number of core parti- 
cles of the core-halo system, we considered an integrated 
MF density profile, f{x) = Airy'^ p(y)dy (see Fig. 0J. 
As it follows from the Figure, the MF core-halo state in- 
deed contains a distinct core with a sharp boundary of 
the radius Xc ~ 10~^ relatively insensitive to the energy 
in the range we considered, |e| < 0.5. Using this MF core 
radius, we located cores in the MD core-halo systems 
which contained very similar to the MF cores number 



TABLE II: Equilibrium MD and MF results for a core-halo 
state at e = -0.339, N = 250, and < r < 1500 



MD 



MF 



e 
a 

..2\2 , 



core 



-0.3392 ±2 X 10" 
2.9 ±0.1 
-1.5 ±0.1 
0.99 ± 0.03 

48 ±2 



-0.339 
2.94 
-1.46 
1 

47.6 
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FIG. 4: Integrated MF density profiles fu{x) of a uniform 
state (dashed line) and fc-h{x) of a core-halo state (solid line) 
for e = ecoii 

of particles (see Table Hljl . Using smaller core radius re- 
sulted in significant reduction in the number of observed 
core particles. A reasonably small over-estimation of the 
core radius did not affect the results of the MD measure- 
ments: we observed that even in the sphere twice the 
core radius the number of particles is only marginally (at 
most by 8%) larger than in the core. 

To check if the system has more than one core, we 
performed search for the second-largest core of the same 
radius Xc- We looked for a largest group of particles 
which are within Xc from a single particle with none of 
these particles belonging to the first, largest core. We 
never observed the second-largest core containing more 
than 2 particles; most of the time it contained only a 
single one. 

In Fig.|Slwe present the MD velocity distribution func- 
tions W{u) for a core-halo and uniform states; shown 
W{u) confirm the MF prediction for the Maxwellian form 
of these distributions. 

As we mentioned in the previous section, we were un- 
able to resolve the high-density part of the radial density 
profile due to the core motion. However, an indirect com- 
parison between the radial distribution of particles in the 
MF and MD was made using the radial distribution func- 
tion. The MF radial distribution function Cmf{x) was 
computed as 

Cmf{x) = -i^ / p(x')p(x + x')dx'. (5) 

The good agreement between the MF and the MD radial 
distribution functions is illustrated in Fig. El This indi- 
cates that the mutual distribution of particles is correctly 
predicted by the MF theory. 
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FIG. 5: The MD velocity distribution functions W{u) of a 
core-halo state with e = —0.339 (solid line) and a uniform 
state with e = —0.3 (dashed line). In both cases A*' = 250. 
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FIG. 6; MF (dashed line) and MD (solid line) radial distri- 
bution functions C{x) of a core-halo state with e = 0.25. The 
step at a; = 1 in the MF C{x) is caused by the localization 
of the core exactly at i = and a sharp boundary of the 
container. 

To summarize, for all the quantities considered, we ob- 
served no systematic deviations between the MF theory 
and the MD data. 



COLLAPSE 

According to the MF theory, if the energy of the uni- 
form state becomes lower than tcoii ~ —0.339, the system 
should undergo a collapse to a core-halo state. To study 
the collapse, we considered several uniform systems with 
the energies ranging between e = —0.5 and e = —0.3. The 
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FIG. 7: Time dependence of the kinetic energy ekin (top) and 
the virial variable a (bottom) of the collapsing uniform state 
with e = —0.5 and A'^ — 125. The dashed horizontal lines indi- 
cate the equilibrium values of ekin and a of the corresponding 
core-halo state. The data is averaged over St — 100 time 
intervals. 
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FIG. 8: Collapse in systems with e — —0.5 and different 
numbers of particles, A*' = 125 (dashed line) and A'^ = 250 
(solid line), shown in relaxation time units, rln N/N [Till . 
Horizontal line shows the kinetic energy of the target core- 
halo state. The data is averaged over St — 100 time intervals. 



systems were initialized according to the MF density dis- 
tributions. For systems with e < ecoii the particles were 
distributed according to the MF density profile for ecoU ■ 
In a perfect agreement with the MF theory, a uniform 
state with e < ecoU undergoes a gradual transition to a 
core-halo state with a typical timescale of TcoU ~ 10'* for 
= 125 - 250 particles. An example of the time depen- 
dence of the kinetic energy and the virial variable for a 
collapsing system is shown in Fig. [7\ We observe that if 
the number of particles is increased but the rescaled en- 
ergy e is kept fixed, it takes generally longer time for the 
collapse to be complete. Our results (see Fig. (H)) quali- 
tatively confirm that the characteristic time for the full 
collapse scales as ^ quantitative study of the 

dependence of the collapse dynamics on the number of 
particles requires much faster simulation code, however. 





10000 



20000 



In the above examples, the energy was set to e = —0.5 
which is well below ecoU ~ —0.339, and as a consequence 
the collapse started immediately at r = in all simu- 
lation runs. If the system energy is Scoii, the noticeable 
increase in kinetic energy and decrease of the virial vari- 
able, characteristic for collapse, start not exactly at t = 
but with a small delay (Fig. ^ which varies from run to 
run from almost zero to about t « 1500. This indicates 
that the MD system is able to overcome the metastability 
at or near ecoii ■ The observed uncertainty is likely due to 
the relatively small number of particles. 

As we increase the energy above ecoU, the stability of 
the uniform state increases which results in a longer life- 
time of such state with respect to collapse. In Fig.^1 an 



FIG. 9: Plots of the kinetic energy ekin (top) and the virial 
variable a (bottom) vs time r for a system with e — ecoii ~ 
—0.339 and A'' = 250. The data is averaged over St — 100 time 
intervals. 



evolution of a system with e — —0.3 is shown. The system 
stays in the uniform state for about 6t « 5000 before the 
collapse starts, after which the evolution proceeds qual- 
itatively similar to the collapses in systems with lower 
energies. 

To compare the temporal evolution of the kinetic en- 
ergy, virial variable, and the number of core particles, the 
relative variables e'fe„(r), <„(t), and A^^ore(^)> all de- 
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FIG. 10: Plots ol the kinetic energy etin (top) and virial 
variable a (bottom) vs time r for a system with e = —0.3 and 
A'^ = 250. The data is averaged over 5t = 100 time intervals. 
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FIG. 11: Plots of the relative values of (from top to the 
bottom) number of core particles A^core (''"), virial variable 
cr^i„(r), and kinetic energy e'^inir) for the system with e — 
-0.339 and iV = 250. 

fined as e^^„ (i) = [efem(i)-efem(u)]/[efcm(c-/i)-efcm(w)], 
are plotted in Fig^] The values €kin{u) and ekin{c — h) 
correspond to the uniform and core-halo states in equi- 
librium. 

FigurelTTlindicates that during the initial stages of col- 
lapse the core grows faster than the kinetic energy and 
the virial variable. In addition, one can notice large re- 
versible fluctuations in the number of core particles (the 
core grows up to 12% of its final value and then disap- 
pears) that are not matched by comparable scale fluctu- 
ations in the kinetic energy or virial variable. All these 
observations suggest that the density evolution causing 
the core formation plays the leading role in the process 



of collapse while the relaxation of kinetic energy follows. 
Once the collapse has started, the core grows to about a 
half of its final size in only Stcotc ~ 10'^ for systems with 
N = 125 - 500 particles, while the changes in kinetic 
energy during this interval of time are small. After this 
rapid initial stage the system relaxes more slowly, and 
finally after TcoU ~ 10^ reaches the equilibrium core-halo 
state. Our observations strongly suggest that the growth 
of the core takes place through a sequential absorption of 
single particles rather than through hierarchical merging 
of smaller cores: We never detected other cores contain- 
ing more than two particles. Although the kinetic energy 
relaxation trails behind the the core formation, the veloc- 
ity distribution function remains Maxwellian throughout 
the whole evolution with the temperature corresponding 
to the corresponding value of the kinetic energy. This 
is caused by the fast velocity relaxation {r-uei < 1) as 
discussed in the previous section. 

EXPLOSION 

It is natural to assume that if a system exhibits a col- 
lapse, it should also exhibit an explosion which is the 
reverse to the collapse transition. According to the MF 
theory, such explosion should take place when the core- 
halo state becomes unstable, i.e., when e > eexpi ~ 0.267. 
To check this prediction, we initialized the MD system 
according to the MF equilibrium core-halo state and fol- 
lowed its evolution. As in the study of the collapse, for 
initial states with e > eexpi we used the MF density pro- 
files of the highest energy locally stable state, i.e. of the 
state with e = ^expi- 

We observe that a system with sufficiently high energy, 
such as e = 0.5 in Fig. El or e = 0.4 in Fig. El indeed 
undergoes an explosion which brings it to the uniform 
equilibrium state. During such an explosion, the state 
variables such as kinetic energy and virial variable con- 
tinuously change from their equilibrium core-halo state 
values to the uniform state ones, and the core gradually 
sheds particles until only one particle is left. The main 
features of an explosion (Figs. IT^ and IT^ resemble those 
of a time-reversed collapse. The kinetic energy evolves 
relatively uniformly, while the number of core particles 
changes only slightly during the first stages of evolution 
and rapidly decreases at the final stages. In the example 
presented in Fig. El the explosion is complete after the 
time texpi ~ 15000, which is noticeably less than the time 
for a collapse tcoU ~ 10^ (see Fig. (HI for a system having 
the same number of particles {N = 250). However, the 
latter is rather vaguely defined due to larger fluctuations 
in a core-halo than in a uniform state. 

Similarly to a collapse, the system remains thermal- 
ized in the velocity space during an explosion. The ve- 
locity distribution remains Maxwellian throughout the 
evolution with the temperature corresponding to the cur- 
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FIG. 12: Plots of the kinetic energy e'fei„(T) (top) and rel- 
ative number of core particles N^orei''') (bottom) (defined as 
in Fig. IIH vs time r for a system with e = 0.5 and A'^ — 250. 
The kinetic energy is averaged over St — 100 time intervals. 
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FIG. 13: Same as in Fig.HSlbut for e = 0.4. 

rent value of kinetic energy. As an illustration, Fig. FT"!] 
shows the ratio of the moments of velocity distribution, 
19{v'^)^ /5{v'^), which should be 1 for a Gaussian distri- 
bution. 

However, as it is evident from a comparison between 
the Figs. [T^ and ICT as e gets closer to eexpi the explosion 
takes longer to initiate. We have observed, that even for 
e = 0.3, which is noticeably larger than eexpi ~ 0.267, 
the explosion does not happen during the first r = 30000 
of evolution. This suggests that either the MF value for 
^expi is incorrect, or during the incitation of the system 
we somehow prepare the system not exactly in the equi- 
librium (metastable) core-halo state. If the latter is the 
deviation from the equilibrium most probably 
takes place in the core, as because of its compactness. 



1.3 



1.2 - 




0.8 I ' ' ' ' ' 1 

5000 10000 15000 

T 

FIG. 14: Plot of he ratio of the moments of velocity distribu- 
tion, 19(u^)^/5{«^), vs. time r for a system with e = 0.5 and 
N = 250. 

its equilibration with the rest of the system may take a 
rather long time. Using the current MD setup, we were 
unable to determine a reason for this apparent discrep- 
ancy. 

CONCLUSION 

In the previous sections we have presented the follow- 
ing molecular dynamics results for the self- attracting sys- 
tems with soft Coulomb potential: 

1. A collapse from a uniform to a core-halo state was 
observed. The timescale for the collapse in sys- 
tems consisting of 125 - 500 particles is of order of 
10^ crossing times and is by the same factor longer 
than the timescale of the velocity relaxation. The 
collapse starts with a fast growth of a core via 
absorption of single particles and continues with 
more gradual relaxation towards an equilibrium 
core-halo state. Metastable states exhibit a finite 
lifetime before collapsing. 

2. A reverse to collapse, i.e., an explosion transition 
from a core-halo to a 

uniform state was observed. The explosion time is 
considerably shorter than the collapse time, being 
of the order of lO'* crossing times (125 - 500 parti- 
cles). An explosion resembles a time-reversed col- 
lapse; the core decrease, which happens by shed- 
ding individual particles, is trailing the kinetic en- 
ergy evolution till the last stages, when the core 
rapidly disappears. 

3. Such molecular dynamics characteristics of the 
equilibrium or the metastable uniform and core- 
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halo states as kinetic energy, wall pressure, num- 
ber of core particles, particle-particle radial distri- 
bution function and velocity distribution function, 
are found to be equal within the statistical uncer- 
tainty of the molecular dynamics measurements to 
the corresponding mean field predictions. 

The long collapse time observed in our simulations 
appears to be an explanation for the apparent discrep- 
ancy between the phase diagram presented in [l^ and 
the mean field phase diagram (see , for example, |lfll|V 
The relaxation time allowed in jl3| before the measure- 
ments of what was considered to be a steady state, 
trei = 3N/\EN\^/^, which is apparently equivalent to 
Trei < 1, is by far insufficient for a system to collapse. 
Therefore, the discontinuities in caloric curves /3 vs e, typ- 
ical for collapse and explosion gravitational transitions, 
were not observed in [13|. 

Although we considered systems only with the soft 
Coulomb potential, we speculate that a likewise simi- 
larity between the mean field and molecular dynamics 
equilibrium properties of the core-halo state exists for 
all "soft" long-range (like a Fourier-truncated Coulomb) 
potentials. This is so because all soft potentials are effec- 
tively longer-ranged than the bare Coulomb one. How- 
ever, the core-halo state in the system with a "harder" 
short-range cutoff may have completely different proper- 
ties from the one considered above, and its mean field 
theory may be inadequate. As for the uniform states, 
their properties are virtually independent on the nature 
of the cutoff (see, for example, [lOj) and their mean- field 
description is universally correct. 

The main goal in the paper was to check the existence 
of collapses and explosions and the validity of the mean 
field data for the self-gravitating systems with short- 
range cutoff. For this goal one or few molecular dy- 
namics runs for each considered system were sufficient. 
However, to be able to study the dynamical features of 
collapses and explosion in more detail and to compare 
the simulations results to various theoretical models, one 
needs to study the relaxation averaged over many ini- 
tial configuration. For example, an interesting question 
is whether a collapse (or an explosion) indeed consists 
of two stages; the first fast stage of collisionless "vio- 
lent relaxation" with particle number-independent rate, 
and the slower second stage of soft coUisional relaxation 
with characteristic time Tr (see, for example, and 
references therein) . Another important question is to re- 
solve the apparent contradiction between the mean field 



prediction for Cexpi and the molecular dynamics observa- 
tions, outlined at the end of the previous section. Such 
studies require a more efficient computation code. The 
main improvement possibly coming from a better force 
calculator that may include various mean-field-like po- 
tential expansions, which are qualitatively justified by 
this study. We leave this for the future research. 
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